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We combine the single site dynamical mean field theory (DMFT) with the non-local GW method. 
This is done fully self-consistently and we apply our formalism to a one-band Hubbard model. Even- 
tually at self-consistency the full self-energy and polarization operator of the system are retrieved. 
Some numerical results, in the metallic as well as the insulator regime, are presented and briefly 
discussed. Depending on the involved interaction (GW) parameters, substantial changes are found 
when the GW self-energy is incorporated. However, the main point of this work is to demonstrate the 
applicability of the method not to make any strict comparison with exact results and experiments. 

I. INTRODUCTION 

The interest for a fundamental understanding of strongly correlated systems has greatly increased, but 
still there is a lack of a satisfactory description. On the other hand, for weakly correlated systems the 
density functional theory^ (DFT) within the local spin-density approximation^ (LSDA), however limited to 
ground-state properties, and the GW approximationMi^ (GWA) suitable for excited state properties, have 
made a substantial contribution to the understanding of sp metals and semiconducturs. Their failure is 
mainly due to a poor description of the strong on-site Coulomb interactions among partially filled d or / 
shell electrons. The insufficiency of the GW method has however encouraged schemes which are all designed 
to treat strong on-site correlations, e.g the LDA-I-U approach proposed by Anisimov and coworkers^ in the 
early 90 's, in order to treat the strong correlations existing in the Mott insulators. There exist several 
similar methodsiSiSiiS that are based on first principles DFT-LSDA Hamiltonians, but the strong Coulomb 
interaction for electrons residing in the localized orbitals are explicitly taken care of via a set of Hubbard like 
parameters, describing static or dynamically self-energy effects. Obviously, there is a necessity to introduce 
in all the LDA+U related methods, a so called double counting correction term for the correlated orbitals2iii. 

Recently, the dynamical mean-field theorj^Mi? (DMFT) has been found to be very successful in the 
treatment of strongly correlated electronic systems. It is a nonperturbative method and has been used 
intensively for various physical properties^"*, such as the famous paramagnetic (PM) Mott- Hubbard metal- 
insulator transition in transition metals, superconducting cuprates, fuUerene compounds as well as organic 
conductors. The DMFT method becomes exact in the limit of infinite spatial dimensionality, and maps 
the original lattice problem onto an interacting dynamical impurity problem, which must be solved self- 
consistently due to its implicit coupling to the surrounding lattice. In the single site DMFT, there is a 
shortage of momentum-dependent or short-range correlations, implying a purely local (on-site) self-energy. 
In the context of spatial ordering and spectral properties that vary across the Brillouin zone, non-local effects 
would of course be crucial. Significant efforts have been made to extend the single site DMFT, to the case 
where the self-energ}ii^iiSiiLi2ii2422iSL22i22i24 exhibits finite-range interactions. The single impurity model is 
replaced by a cluster-impuritjJ^', giving rise to short-range correlations ranging to the boundary of the given 
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cluster. The general idea is that the cluster captures, albeit the finite correlation length, the correlations 
within the original infinite lattice. Some of the approaches, however, breaks the translationally invariant 
nature of the original problem, a scenario not present in the single site DMFT. The corresponding impurity 
problem is considerable harder to solve with the increased number of local degrees of freedom. Present 
techniques are based on the non-crossing approximation2S*2& (NCA), the iterative perturbation theorjii^ 
(IPT) , the Quantum Monte Carlo method^^ or exact diagonalizationiSi2&. An interpolative approachpSS has 
recently been suggested, where a simple pole expansion of the self-energy is used and the unknown parameters 
entering is determined using a chosen set of constraints. 

More recent and probably one of the most promising first principles scheme is the so called "LDA-I-DMFT" 
approacbiSi^Si^, despite the fact that the interaction term for the localized electrons still has to be 
parametrized and the double counting term remains. The parameters are, at least in principle, obtain- 
able from an independent calculation such as e.g the constrained LDA mcthod22iSi24 or from experimental 
data. Note that the screening in the system is not determined from first principles. The feasibility of 
the approach has indeed been demonstrated in the pioneering work by Savrasov and coworkers in the case 
Plutonimn (Pu)'^^ and more recently in a number of other cases^SiSi. 

It is generally believed that the GWA quite adequately describes the long-range part of the screening. 
Short-range correlations, on the other hand, is not taken into account properly by the random phase 
approximation^^ (RPA), however captured by the DMFT approach. Contrary to DMFT, the GWA is a 
perturbative method. The self-energy is given by E = GW, where W is the screened Coulomb interac- 
tion and G is the full Green's fmiction. The frequently used RPA screening {W — Wq) and the zeroth 
order Green's function (G = Go) provide quasi-particle spectra of most semiconductors and insulators 
as well as bandgaps in good agreement with experiment^. However, there is the important issue of self- 
consistenc5*22ii2iiiii242iiiiiSiiSi^. If the GWA should be conserving—, the self-energy requires the Green's 
function as well as the the screened interaction to be evaluated self-consistently. 

The aim of this paper is to combine, fully self-consistently, the GW method with the single site DMFT, and 
present numerical results for a one-band Hubbard model. The "DMFT+GW" approach, recently proposed 
by Biermann et al^, includes no Hubbard-like (parameter) interaction and consequently there is no need 
for the ambiguous double counting term. The main idea is that the large on-site part of the self-energy is 
calculated using DMFT and the off-site (long-range) contribution is taken from the GWA. We will present 
results using a various degrees of self-consistency for the GW self-energy. Related work along this line can 
be found in Refsi^S*^. 

We will study two sites per unit cell in one (ID chain) and two dimensions (2D plane), in order to be 
able to study the formation and stability of different magnetic structure. We solve the single site impurity 
problem using the exact diagonalization method^. In addition to the impurity self-energy, a two-particle 
correlation^^ function is calculated, needed for the evaluation of the impurity screened interaction. Thus, 
the iterative loop will include two quantities to be determined self-consistently: the bath Green's function Q 
as well as the bath effective interaction U. We like to stress that the effective Hubbard U is not a parameter, 
it is in fact found self-consistently. 
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In Sec. II. we describe the method for the calculations. In Sec. Ill we present and discuss the results, and 
in Sec. IV we give a short summary. 

II. THEORY 
A. Single site DMFT 

In this section we establish the necessary concepts and formalisms for the so called single site DMFT, 
a scheme that later on is combined with the GWA. We will consider the Hubbard model with an on-site 
interaction U and nearest neighbor hopping t (t = 1 while leaving the on-site interaction U variable). The 
unit cell will contain two sites, but a generalization is straightforward^^. The model and the corresponding 
Green's function reads 

H = -t al^,^^anja + uYn^i^nraii (1) 

7ni,nj,(T mi 

G™,„,,.(t) = -0(r)(a„„.(r)a,t^.,(O)) + 0(-r)(ai^.,(O)a™.(r)). (2) 

We define positions in the lattice by Rmi = T™ -f- r^, where labels sites within the unit cell and m a 
particular unit cell^^. Using the equation of motion for G (with operator K — H — fJ,N) and assuming a 
local self-energy, T,mi,nj,a — ^ia-SijSmn, one can show that 

G,^-;^(k;w„) = (w„ + + /i^ (k) - Sio.(w„)% (3) 

where the kinetic energy matrix is given byS4 

2{cosk^ + cosky)\ 

~ \2{cosk.^ + cosky) J- ^' 

We have also defined the real space transforms as 

n m 

G^^,nJA^^l) - ^ ^e''"^""G,,a(k;ztuz)e-^'^-'^- (6) 

k 

where the lattice has N unit cells. In the Matsubara formulation, wc adopt the definition 

r/5 



Jo 



(7) 

n 

where z/„ denotes the Matsubara (odd) frequency for fcrmion propagators. For bosons we use Un (even) as 
a convention. 

= (9) 



Inversion of Eq. Q gives the lattice Green's function 

. X /w„ + M - S2<t(w«) -2{coskx + cosky)\ 

where I?(k, = + — S2cr(w„))(ii^„ + — Sicr(it'ri)) ~ ^(cos/cj; + cosA:y)^. The local (impurity) Green's 
function is calculated using the diagonal elements; 

Gicr(ii'n) = -^^Gii<T(k;w„). (12) 

k 

Regarding the corresponding self-energy, we remark that in the case of single site DMFT, no causality 
problems occurs, the lattice self-energy is identical to the impurity self-energy: Ytij^Ck^iVn) — '^ia{ii^n)Sij. 
At DMFT self-consistency the Green's function calculated using Eqs. H11I12|I must coincide with the one 
extracted from the impurity model. We have determined the site and spin dependent impurity Green's 
function using the exact diagonalizationS^ (ED) Lanczos method for the single impurity Anderson model. 
In the present case (zero temperature; /3 oo), we have solved an effective impurity model for each site 
i = 1,2, given by 

<r k=l k=l 

+ Un^^nil (13) 

where Sd = —fJ- is the energy of the localized level on the impurity site. The second term gives the energy of 
all the bath (conduction band) electrons, which are labelled by k — 1, ..,Ns — 1. The hopping between the 
bath states and the impurity state is described by the third term, where Vika is a hopping matrix element. 

The DMFT approach maps the original lattice problem defined by the Hubbard model onto a self-consistent 
solution of the Dyson equation in Eq. and the (auxiliary) impurity problem defined by the bath Green's 
function 

gr/ (j^„) = G,r; {ii^n) + S,. (*i.„) . (14) 

In order to initialize the iterations it is sufficient to guess the parameters of the Anderson model, Sika and 
Vikcr, as well as the bath Green's function. We construct 

Qtcy{Wn) = ^[(it'n + M)l5y + /ly (k) - -Bjcr(w«)(5ij]~^ (15) 

k 

where Bi„ is a chosen suitable external field (in the PM case Bi^ — 0) . Solving the effective impurity model 
we derive the self-energy Sig. = Q^J^ — G^J and proceed with the inversion of the matrix in Eq. ||2J). Finally 
we update the bath Green's function using Qia{ivn) — [1/ X^k Guai)^', ivn)lN -t- Y^ia{wn)]^^ and mix it with 
the previous one. The bath Q^^ is represented by the U — Q impurity Green's function: 

Ns-l y2 

k—1 
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in order to provide us with a new set Anderson parameters, found by a fitting procedure. The best choice 
is found by minimizing the functionSSiS 

= E \Q^.\Ns,^yn) - Ql.\^yn)\|y^ (17) 

n— 

for each site i and spin channel a. The convergency with respect to N s is very fast. We found that already 
A'^s — 1 = 7 bath states is sufhcient to describe the continuum of conduction states. The DMFT cycle is now 
closed: at hand we have a new set of Anderson parameters (which defines the impurity problem) as well as 
an updated bath Q . At self-consistency, the Green's function from the impurity problem should be equal to 
one obtained from summing the momentum-dependent lattice Green's function over the Brillouin zone, as 
done in Eq. l(T^ . 

When DMFT is combined with the GWA, the the impurity charge response is entering the formalism. 
The two-particle response is defined by 

= -{p,{T)p.,)e{T)~{p,Mrm-T) (18) 

where Piir) = fii{T) — rii, Pi{T) = e^^'^ piC^^^'^ and the total charge on the impurity is denoted by rii = 
ni] + nil- From a numerically point of view, the charge response is evaluated on the same footing as the 
Green's function, with the aid of Lanzcos algorithm. All calculations are done for a fix chemical potential 
/i. The total number of electrons in the cell 

n = i ^ n,^ (19) 
is then allow to adjust self-consistently. In the PM case (no doping) ni^ = 1/2 for all sites and spin-channels. 



B. DMFT combined with the GWA 



We now consider a schemeAS, that properly adds the momentum-dependent GW self-energy to the local 
DMFT self-energy, giving rise to a lattice self-energy which describes, in addition to local effects, also long- 
range correlations. The RPA will be used for the screened interaction, implying = G(l — UP'~^^' )^^U , 
where we used ^(r — r') = U6{r — r') for the bare Coulomb interaction. Note that even if v is short-ranged, 
W can have off-site components coming from p'^^ , 

The polarization operator (bubble) in the GWA is given by 

/^«>^(q;zc.,„) = ^/^fr(q;zc.„,) (20) 

where 

i^^^(q;jw^) = ^^^^Gija{q_ + 'k;iiVm + iVn)Gji^{\i\iVn). (21) 

The sum over k includes N = Nk points in the first Brillouin zone (BZ), and q belongs to the irreducible 
BZ. The Green's function in Eq. (|21|l is obtained by inverting the matrix 

Gj^-^(k; iun) = {ivn + p)5ii + hijO^) - Syo.(k; w„) (22) 



where E,yo.(k;w„) is the proper lattice self-energy (see Eq. (|25|l ). In the first iteration, however, the local 
impurity self-energy Eicr(i^'ri) is used. 

The screened interaction fulfills W ^ v + vXv = e~^v, which in the case v{r — r') = U5{r — r'), using 
e = 1 — vP^^ transforms to 

Wij{q;iu;„i) = C/IIy (q; zcj„i) (23) 

where 11 is the matrix obtained by inverting the dielectric matrix [Sij — UP^'^ {q;iijj„)]. If the Coulomb 
interaction v takes into account nearest {V) and next nearest neighbors interaction {V') the dielectric function 
and screened interaction reads in ID and 2D case respectively 

_ fl- UPg^ - 2V'cos{2q)Pg^ -UPg^ - 2Vcos{q)Pg^ 
^ ~ ^UPg^ -2Vcos{q)Pg,^ l-UPgg^ ~2V'cos{2q)Pgg^ 

W = ( + 2V'cos{2q)e^' Ue'' + 2Vcos{q)e^,' 
{ Ue^,^ + 2Vcos{q)e^^^ Ue^g +2V'cos{2q)e^^ 

1 - UPg,"^ - W'cos{q,)cos{qy)Pg^'^ -UPg"^ - 2V{cos{q,) + cos{qy))Pg^'^ \ 
-UPgg^ - 2V{cos{q,) + cos{qy))Pg^'^ 1 - UPgg^ - AV cos{q,)cos{qy)Pg^'^ ) 

_ , Ue^l +Wcos{q^)cos{qy)e'{l Ue^^ + 2V{cos{q:c) + cos{qy))e^l\ , , 



Ut2i + '^V {cos{qx) + cos{qy))e22 ^^^22 + 4T^'cos(gx)cos((7y' 
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Like the polarization bubble, the screened interaction is a real valued function on the imaginary axis (even 
Matsubara frequencies) and the diagonal part (i = j) is positive and approaches the bare U for large tJm, 
implying that the correlated part (frequency dependent) of W goes to zero {Wfgq;,iu}m) ~ Sij/iiUmY when 
iOm 00 and V — V — 0). 

Finally we achieve for the GW self-energjiS (5]j^^(q; w„) = C/ni_cr5ij-l-S-j^(q; ij^„) and W^^ = Wij-USij) 

Sy,,(q;ji^n) = -^^-^^G.y,,(q-k;iz/„ -itj„,)Wfj.(k;itj„,) 

= ^]]^^Y1 Gy<x(q-7^k;^^^/„-^tJ™)Wi;;■(k;iw™) (25) 

where W{TZ\<i) = VF(k) has been used. TZ denotes a rotation matrix corresponding to a point-symmetry 
operation^. The particle number used for the Hartree-Fock part (S^-'^ = Uni^cr^ij) is calculated using 
the impurity Green's function, however at self-consistency the impurity Green's function and the local one 
should be identical (the k dependent lattice Green's function summed over k). The total lattice self-energy, 
corrected for double counting, and to be used in the construction of the next G^^ can thus be written as 

Ys,j„{q:,ivn) = Y.f^^{q:,iVn)-dtj^^Yjf^{\i;iVn) + '^ia{iVn) 



k 



= S«^(q;z^.„)-5,, Y.f^{^^■^u^)w^ + Y.,„{^v^) (26) 

kelBZ 

where Wk is the weight of k in the IBZ. Note that the local part of S (Sio-) is usually much larger in 
magnitude than the non-local part given by [E*^^ — l/iVfc X^k ^'^^]- 

Finally the local G to be used to find the bath Q via the self-consistency relation: 

{ivn) = {iv^) + E,, (^f.„) (27) 
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can be written as 

Gicr(ii^„) = — ^ G.ua(k; w„) = ^ GMo-(k;w„)wk (28) 

k kGlBZ 

where the diagonal-elements Gna{k.]ii'n) are found from inverting 

G~-^(k; ivn) = {iun + iA^ij + %(k) - Sya(k; iun) (29) 

with the self-energy from Eq. I|26|l . 

In an ordinary single site DMFT calculation the impurity problem is solved for fixed on site U and only 
the bath Q is updated and determined self-consistency via Eq. 127|l . It is however desirable to solve the 
impurity problem with an updated or an effective Hubbard interaction. The static impurity charge response, 
Xi(ict;„i = 0), is used to construct the static impurity screened interaction and polarization 

W^iiuJr,^ = 0) = U,+U^Xt{ii^m = 0)U, (30) 

P,{iiOm ^0) ^ Ur^~ Wr\iu;m = 0) (31) 

where Ut is the effective Hubbard onsite-interaction used for the solution of the impurity problem at site i. 
Then the full polarization kernel can be written, using Eq. H21|l . 

P^J{q;^UJ^) - i^^'^(q;iw„)-5,,^^i'^^(k;iw™) + P,(iu;,„ = 0) 

k 

= P^^(q; ILUm) - % P?^^^'^ i^m)wk + P^U^rr. = 0) (32) 

kelBZ 

a relation analogous to Eq. H2t)|) . Then the local screened interaction reads 

-^"YWu{k;ioj,n^O)= Wu{k;iuJm = 0)ww (33) 

k kelBZ 

where the diagonal-elements Wa (k; icurn = 0) are found from inverting 

W,-i(k;ic^™ = 0) = U-'6,j-P,,{k;iiu^ = 0). (34) 

Note that here the bare U is used (same for all sites Ui ~ U). In the case when hopping to neighbors is 
allowed we have to substitute the diagonal term U^^Sij with the inverse of the bare Coulomb matrix 

( U + AV'coskxCosky 2V{coskx + cosky) \ ,„_^ 
y 2V{cosk.j; + cosky) U + AV' coskxcosky J ' ^ ' 

Finally the updated effective interaction is found from the self-consistency relation 

= l/Y. ^"(k; = 0)/A^fe + m^;,^ = 0) (36) 

k 

a relation analogous to Eq. H27|) . however, only the static value oilA^^ is used in the solution of the impurity 
problem. 

The spectral function is given by 

A(k; c.) = - ^Im ^ ^ G.,. (k; a;) (37) 



where Nt — 2 in the antiferromagnetic (AF) case. In order to obtain the Green's function on the real energy 
axis, we use the Pade approximation for the self-energy Sij^Ck; w„) Tiij„{\<i]Lj). 

In order to study the stability of different phases, the total energy (per site) is calculated using 

E = Tr{/^(k)G(k; lUn)] + ^Tr{S(k; iz/„)G(k; lUr,)] (38) 



where 



^Ei^Ei^EE (39) 



The hopping and self-energy matrices are given in Eqs. I|24l26f) and the the Green's function matrix is found 
by inverting Eq. (|29|l . 

C. Computational details 

Some care has to be taken when performing the Matsubara sums for the polarization bubble and the 
self-energy in Eqs. H21I25() . The bubble can be written as 



Pi^j^ {^.\i^ni) = ^ ^ ^^[G.y>(q + k;iw,„ -I- iz^„)Gji„(k; w„) -I- 

GiJa{^. + '^■■^'^rn + iV-n-l)G*i„{\fi;Wn)] (40) 



n>0 k 



where we have used that G{iv-n) = G* {ivn-if'^ ■ The polarization is real valued on the imaginary axis 
(even Matsubara frequencies) and the diagonal part (i = j) is negative; Pj*^'^(q; ioJm) ^ Sij/{iujm)'^ for large 
LOm- We note that for large n the first term behaves as 

(3 E j(t^ + Pn)iVn ' ^^"^^ 

n ^ ' 

This is however not the case for the second term, but yet we find that the following procedure is appropriate: 
If the Matsubara sum is done for all frequencies on the imaginary axis the result is ^ji j ^ for m = 0, otherwise 
zero. We have subtracted the term ^ijl fi'Ylm l/*('^r?j. + Vj^iv^ in Eq. (|40|) (where, of course, the sum is done 
for finite n) and consequently added —SijP/A. The second term, however is large whenever (m — n — 1) 
around zero, due to G(m — n — 1), even if G*{n) is decaying for large n. Therefore the upper limit for the 
n-sum in Eq. (|4()|l is chosen to depend on m. Thus we evaluate 

P^^i^l'^^'^m) = ^ E Jl-^[GiJ'y(.<l+'^',^^rn+i'^n)Gj,„{k;il'n) + 
^ n>0 ^ k 

Gya(q + k; iujrn + w_„_i)G*,^ (k; ji/„) 

- ^^^h ^ V + -( ^ ^ u- yJ ^ - K^m)S^iP/^- (42) 

where Np{m) = Np + m. The polarization is calculated as described above for m = 0,Nh, whereas for 
m = Nh + l,Ng we fit to 

P^'^ici-^zu;^) = S,,-^ (43) 



where P° is a positive constant chosen for a continuous match. 
The correlated GW self-energy is given by 

Gij„{ci-']lk-,iu„+iu;rn+i){W%{\<i;iu;rn+i)]- (44) 



/3 ^ -/Vfe 

^ m>0 n kGlBZ 



The first term, is large whenever (m — n) aroimd zero, due to G{m — n), even if the screened interaction 
is decaying for large m. This means that the upper limit for the m-sum should depend on n. We have 
performed the sum for m = 0, Ns{n) where Ns{n) — Ns + n. 

The impurity (Anderson Hamiltonian) is solved using the updated effective U, not the bare U. To be 
consistent, the localized level in the impurity model is updated using Ed = —fJ^ = + A/i). In the half- 

filled case A/i = (hole doping A// < 0). We have scaled the bath Q^^ as well as the impurity self-energy 
T,icr. We have 

<^ija(l*;*^n) = {it^n + IJ')Sij +hij{\s.) -Y,ij„{]s:;iUn) 

= [iun + U/2 + A/x - T.i„{ivn))5ij + /ii,(k) - Sg^(k; iu^) 

= {iVn + A/i - [^i^{iVn) - U/2])5ij + hij (k) - Sg^(k; iVn) 

= {ivn + A/i - tia{ivn))kj + /lii(k) - Sgr(k; ivn). (45) 
The GW self-energy includes the double counting term. We also have 

QiJ'{Ns,ivn) = {ivn-Sd)- 

fe= 

Ns- 

= [ivn + m) - 

= (ii/„ + U/2 + A/i) - — . (46) 

fc=l 



Ns-l 2 



Thus 



Qi,\N,,iv^) = (*:.„ + A/i) - ^ . (47) 
with = — U/2. The self-consistency relation: 

gr\ivn)-U/2 = Gr^{iun) + ^ia{iyn)-U/2 

gr'iii^n) = Gri(ii.„) + S,<,(iz.„) (48) 

We note that the Hartre-Fock (impurity) self-energy can be written as 

i:f/=Uni.„. (49) 

In the half-filled case tii-^ = 1/2 for all sites i and spin- channels, so TiialiVn) is the impurity self-energy 
with the static Hatree-Fock part removed. 
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Iterative steps: 

1. For each site i in the unit ceU and spin-channel a, we have to solve an impurity problem. The 
Anderson Hamiltonian, which is defined by sa, {en-cr, Vik^} and the effective Hubbard Hi, is solved in order 
to get the impurity Green's function Gj^. and the static response Xii^'^m = 0)- Using the response we can 
calculate the screened interaction for the impurity; W, {ium = 0) = Wj + UiXi {ii^m = 0)ZY, as well as the 
impurity polarization Pi{iLOm = 0) = Uf^ — W^^{iu)m = 0). 

2. Derive the (scaled ) impurity self-energy from Sjo-(«^'n) = Gi^^iii^n) — G^a{wn)- Here we use the 
bath Green's function from the previous iteration. In the first iteration we have to guess the Anderson 
(bath) parameters as well as the bath Green's function. 

3. With the impurity self-energy we construct (S^^(k;iz/„) from the previous iteration which includes the 
double counting term for i = j) 

G-a(k;zi/„) = {^Un + A^l-ti,{iUr,))SiJ+hiJ{k)-J:^'^'{k■,^lyn). (50) 

Using Gija{\i;ii'n) we construct the updated GW self-energy to be used in the next iteration and then we 
calculate the local Green's function X^j^ G,io-(k; zz/„) using the impurity self-energy and the updated GW 
self-energy. We also calculate the local screened interaction W,i (k; iw^ = 0). 

4. Update bath Green's function using Qiaiii'n) = [l/X^k^«»o-('*^'^'^") + ^^'^ effec- 
tive interaction using Hi = [1/ ^aO^' = 0) -|- Pi{icOrn = 0)]~^. 

5. Mix old (bath Q used in step 2) and new (bath Q from step 4). Same mixing for old effective 
interaction ( U used in step 1) and new ( U from step 4). 

6. The mixed bath Green's function is then fitted {Qi„^{ii') ~ Q^J^{Ns,iv)) in order to deter- 
mine the updated parameters {cika, Vik^}- 

7. Now we have a new set of parameters (which defines the impurity problem) so we go back to 
step 1. We also have a new bath Q to be used in step 2. At self-consistency the Green's function obtained 
from the impurity problem is equal to local one obtained from '}2]^Gn„(yi;ii'n) and the impurity screened 
interaction is identical to ^j^. W,i(k; = 0). 
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III. RESULTS AND DISCUSSION 

We use a simple model system as a test of the feasibility of the method and therefore consider a one-band 
Hubbard model. It is worth to point out that we are mainly interested in how properties, derived using the 
DMFT, changes when GW effects are incorporated as well as the the stability of the iterative procedure. At 
self-consistency we have access to the full self-energy and polarization operator as well as Q and U. In this 
work we focus on the PM solution at half-filling (one electron per site) but not to close to the metal-insulator 
transition. We believe that a careful analysis of the fictitious temperature and the number of bath sites is 
not so crucial when the system is quite far from the metal-insulator transition. All results presented here 
will be for four bath-sites {Ns — 4) . The system studied consists of two sites in the unit cell (denoted 1 and 
2) and we impose no constraints on different sites and spin-channels i.e in the paramagnetic case we will 
obtain four identical solutions. If the system initially is in the metallic PM phase, the system can during the 
iteratively procedure, end up stable in the insulator AF phase when convergency is reached. Such a scenario 
is of course not possible if only one site and spin is considered per unit cell. We will assume that all energies 
are given in eV. 

A. ID chain 

Although a Luttinger liquid we will consider the ID chain (bandwidth 4) and we have chosen U — 2 and 
U ~ lA as prototypes for a metal and an insulator respectively. We have checked the convergency with 
respect to the number of bath-sites. In Fig. ^the imaginary part of the on-site lattice Green's function is 
displayed as a function of imaginary (odd Matsubara) frequencies corresponding to the inverse temperature 
(3. Convergency test with respect to the number of points in the (IBZ) in addition to the energy-range 
parameters {Ng, Nh, Np and Ns defined in paragraph II C) has been performed as well. We will first 
discuss a typical metal. Apart from the on-site interaction (short-ranged) U the present GW approach also 
contains the off-site (long-ranged) interactions V and V' (see Eq. I35|l . Quite naturally the significance of 
the GW effects are in some sense tuned by the magnitude of these off-site interactions. We have chosen the 
parameters V = 1.5 and y = 1.2 in the metal case U — 2. This choice of parameters are not, at this point, 
dictated by any physical grounds. However we believe that parameters chosen are in a such a range that at 
least some comprehensive statements can be made. The difference between using /? = 10 or /3 = 20 is very 
small (the number of Matsubara energy points in the low temperature case was increased correspondingly) 
and if not stated otherwise the inverse temperature is /? = 10. 

In Figs. H2I3|I the /c-dependent lattice Green's functions are shown in the low-energy region. In the DMFT 
case the total self-energy is merely composed of the local impurity (fc-independent) self-energy defined in Eq. 
14 (S*^^ = 0). Obviously the inclusion of the GW self-energy is quite substantial for small energies. We like 
to stress that the total self-energy (Eq. 126(1 exhibits non-diagonal site contributions originating from the GW 
kernel, influencing the Green's function and consequently the spectral properties. The displayed behaviour 
of the Green's function has been observed by several authorsSSiSiSiiS. Capone et al.^ have found, in the 
metallic region, that the inclusion of a cluster DMFT approach will give rise to a dip in the in the imaginary 
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part of the on-site Green's function, albeit characteristic of an insulator. 
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FIG. 1: Imaginary part of the site-diagonal Green's function at the F-point for ATj = 4 and 6. Parameters used for 
DMFT-I-GW: Nk = 33, Ng = 512, Nh = 128 and Np = Ns = 64. 
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FIG. 2: Imaginary part of the site-diagonal Green's function at the F-point. Parameters used for DMFT-I-GW: 
Nk = 33, Ng = 512, Nh = 128 and Np = Ns = 64. 
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FIG. 3: Real part of the site-nondiagonal Green's function at the F-point. Parameters used for DMFT+GW: 
Nk = 33, Ng = 512, Nh = 128 and Np = Ns = 64. 



The GW derived polarization (Eq. I4UII and screened interaction (correlated part) (Eq. I23|l are displayed 
in Figs. H4I6(I as a function of imaginary (even Matsubara) frequencies. Note that if the full polarization in 
Eq. 1321 is required one has to correct for double counting and merely add the static impurity contribution 
{Pi{iuj,n = 0) = -0.47). For large Matsubara energies the diagonal part approaches 2V' and the non-diagonal 
part 2V as can be derived using Eq. H24|) . which is numerically confirmed. In Fig. Q we show the screened 
interaction at the F- point along the real axis using the Fade approximation. We observe that the static value 
ReM^(O) is merely a constant below the main excitation peak and slightly larger in the case of non-diagonal 
screening. However it is well-known that in the RFA the screening is overestimated at short distances. From 
a physical point of view this fact is easily understandable: a positive hole is surrounded or screened by a too 
tightly drawn electron cloud, due to the fact that exchange and correlation effects are neglected among the 
screening electrons. 
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FIG. 4: Real part of the polarization function at the F-point. The static impurity contribution is Pi(iuim = 0) 
-0.47. Parameters used for DMFT+GW: Nk = 33, Ng = 512, A^^ = 128 and Np = Ns = 64. 
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FIG. 5: Real part of the site-diagonal (correlated ) screened interaction at the F-point. The bare Hubbard U has 
been subtracted. The impurity screened interaction is 0.5 and the effective Hubbard Vl = 0.7. Parameters used for 
DMFT-I-GW: Nk = 33, Ng = 512, Nh = 128 and Np = Ns = 64. 
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FIG. 6: Real part of the site-nondiagonal screened interaction at the F-point. Parameters used for DMFT+GW: 
Nk = 33, Ng = 512, Nh = 128 and Np = Ns = 64. 



The self-consistent values of the impurity screened potential and the effective Hubbard on-site interaction 
(both defined in Eq. I30|l were found to be W{iuJm — 0) ^ 0.5 and U — 0.7 respectively. Thus at self- 
consistency, the effective impurity problem offers an on-site interaction which is more than a factor of two 
smaller than the bare U = 2. For illustration we show the charge impurity response function along the real 
axis in Fig. |S1 derived using the effective Hubbard U = 0.7. 
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FIG. 7: Real and imaginary part of the site-diagonal screened interaction at the F-point. We used an artificial 
broadening of 0.5. Parameters used for DMFT-I-GW: Nk = 33, iV,, = 512, Nn = 128 and Np = Ns ^ 64. 
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(0(eV) 

FIG. 8: Real (solid line) and imaginary ( dashed line) part of the impurity response function (site 1) defined in Eq. 
CHlfor the effective impurity. Parameters used for DMFT+GW; Nk = 33, Ng = 5f 2, Nh = 128 and Np = Ns = 64. 



As discussed previously, in the single site DMFT case the solution to the impurity model is extracted 
using the bare Hubbard U — 2, however in the DMFT+GW scenario the impurity model is solved with 
the effective (weaker) interaction lA. The magnitude of the impurity self-energy scales with the size of the 
on-site interaction making it somewhat cumbersome to compare different impurity self-energies obtained 
with different on-site interaction strengths. However the quantity one really should compare is the total 
self-energy entering the theory i.e the U = 2 impurity DMFT single site self-energy should be compared 
with the full self-energy in Eq. H26|l . In this work a critically comparison will not be done, we briefly discuss 
spectral properties, which however strongly depends on the self-energy. 

Prior to the discussion about spectral properties we intend to make some statements about the derived 
self-energies. In all the cases studied we have observed the characteristics of a metal or an insulatoui^: 
lmS{iuj) < increases linearly or diverges when a; ^ 0+ respectively. Regarding the magnitude of the GW 
self-energy it depends strongly on the fc-point, but in general the non-diagonal KeT,f^ is quite large and 
ImE]']^ is smaller (in comparison with relevant quantities). 

It is a delicate matter to extract real frequency dynamical information from imaginary axis data. The 
commonly used quantum Monte Carlo impurity solver uses maximum entropy based methods^^. In the 
present work, adopting the Lanczos routine for solving the impurity problem, we used the the Fade approx- 
imation when performing the analytical continuation. For example in order to obtain the GW self-energy 
and spectral functions we must do an analytical continuation from the imaginary axis {iuj (uj + i5)). 
However, the impurity self-energy on the real axis can be extracted using the self-consistency relation in 
Eq. 03 The corresponding resuhs are shown in Figs. (HmH) (DMFT) and Fig. UHl (DMFT+GW). In the 
DMFT+GW scenario the impurity is solved with an effective Hubbard interaction U = 0.7, reflected in a 
substantial reduction in the magnitude. As a comparison with Fig. the self-energy derived using the 
Fade approximation is displayed in Fig. 1121 
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In the metallic case the self-energy exhibits the Fermi-liquid behaviour: ReE(cj) ^ (1 — 1/Z)uj (or equiv- 
alently lmS{iuj) ^ (1 — 1/Z)uj) and Iml](ijj) ^ — tj^ for uj close to zero where 

^=(l-^^Uo) (51) 

denotes the quasiparticle renormalization factor. In the insulator case the slope of Rel](cij) changes sign 
(ReI](Lj) l/uj for uj ^ 0) and lmT^{uj) is peaked at the chemical potential and zero in the gapSfi as evident 
from Fig. lfTT|l . 
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FIG. 9: Real and imaginary part of the impurity self-energy (site 1 and spin up) for a typical metal. We used an 
artificial broadening of 0.75. Parameters used for DMFT: Nk ~ 33 and Ng = 512. 




FIG. 10: Real and imaginary part of the impurity self-energy (site 1 and spin up) for a typical metal. We used an 
artificial broadening of 0.75. Parameters used for DMFT-I-GW: Nk. = 33, Ng = 512, iV;, = 128 and Np = Ns ^ 64. 
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FIG. 11: Real and imaginary part of the impurity self-energy (site 1 and spin up) for a typical insulator. We used 
an artificial broadening of 0.75. Parameters used for DMFT: A^'fc = 33 and Ng = 512. 



DMFT impurity self -energy 


10 Pade coefficients 




Broadening . 05 


eV 






















U=2 

; p =10 

i N, = 4 



-20 20 

co(eV) 



FIG. 12: Real and imaginary part of the impurity self-energy (site 1 and spin up) for a typical metal. Parameters 
used for DMFT: Nk = 33 and Ng = 512. 



We next discuss spectral properties. In order to achieve the local density of states (LDOS) we have solved 
the impurity model on the real axis and then extracted lmGia{i^)- As can be seen in Fig. E| the LDOS 
is symmetric (half-filhng n = I) and shows the typical Fermi metallic characteristics; a quasiparticle peak 
surrounding the two Hubbard bandsi^. 
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FIG. 13: Local density of states for site 1 and spin up. The chemical potential corresponds to energy zero. We 
introduced an artificial broadening of 0.25. Parameters used for DMFT+GW: Nk = 33, Ng = 512, A''^ = 128 and 
Np =Ns = 64. 



With aid of the Pade approximation and Eq. 03 we calculate the spectral functions. The zone-center 
spectral function is visualized in Fig. (|14(l . A significant change of the quasiparticle peak position is clearly 
seen at the F-point, where the downward shift is around 0.4. The corresponding dispersion in the T-X 
direction is displayed in Fig. 1161 Interestingly when only one iteration with the GW kernel is performed on 
top of a self-consistent DMFT calculation, the dispersion essentially coincide with the DMFT one. 
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FIG. 14: Spectral function at the F-point. The chemical potential corresponds to energy zero. The solid line 
corresponds to the DMFT-I-GW case and the dashed line to the DMFT case. We used an artificial broadening of 0.5. 
Parameters used: A^^ = 33, Ng = 512, Nh = 128 and Np = Ns ^ 64. 
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FIG. 15: Spectral function at the X-point. The chemical potential corresponds to energy zero. The solid line 
corresponds to the DMFT+GW case and the dashed line to the DMFT case. We used an artificial broadening of 0.5. 
Parameters used: Nk = 33, Ng = 512, Nn = 128 and Np = Ns ^ 64. 



> 




-2 t ♦ 



V'=1.2 
V=l .5 

U=2 

♦ ♦ 

□ 



♦ 

♦ □ 



□ 

♦ DMFT 

□ DMFT+GW 

» DMFT + 1 iteration GW 



FIG. 16: Quasiparticle dispersion in the F-X direction. Parameters used: A'^fe = 33, Ng — 512, A*';i = 128 and 
Np =Ns = 64. 



To obtain a notable effect with the GW kernel, one has in general to include the long-range part of the 
bare Coulomb potential and consider nearest {V) and next nearest neighbors interaction (V). For example 
the parameter-set V = V = gives a quasiparticle peak position shifted only by 0.1 compared to the DMFT 
situation (the shift is 0.4 with V ~ 1.5, V = 1.2) which is realized from Fig. El If one compare the lattice 
Green's function in Fig. Hand Fig. UHl it is obvious that DMFT and DMFT+GW with the long-range part 
excluded {V = V = 0) are quite similar. 
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FIG. 17: Spectral lunction (DMFT+GW) at the F-point . The chemical potential corresponds to energy zero. The 
soHd hne corresponds to V — 1.5, V' = 1.2 and the dashed hne to V — V' = 0. We used an artificial broadening of 
0.5. Parameters used: Nk = 33, Ng = 512, iV^ = 128 and Np = Ns ^ 64. 
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FIG. 18: Imaginary part of the site-diagonal Green's function (DMFT+GW) at the F-point. Parameters used for 
DMFT+GW: Nk = 33, iV,, = 512, Nn = 128 and Np = Ns ^ 64. 



Next we briefly consider a ID insulator using the parameter-set U = lA, V = 3 and V = 2. In contrast 
to the metal case the screening is less effective giving the self-consistent impurity screened interaction to 
be 13.7 and the effective Hubbard U — 13.9. The similarities between the screened and bare interactions 
indicate that the off-site hopping parameters V and V are too small to give rise to a notable effect, which 
is indeed confirmed by the spectral function shown in Fig. 1191 Furthermore it is worth to note that in the 
strong insulator case the imaginary part of the (site-diagonal) impurity self-energy is diverging for small 
ioj (I](iw) making at least the significance of the diagonal GW self-energy negligible. However 
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non-diagonal GW contributions can influence the spectral functions. 




FIG. 19: Spectral function at the F-point. The chemical potential corresponds to energy zero. The solid line 
corresponds to the DMFT+GW case and the dashed line to the DMFT case. We used an artificial broadening of 0.5. 
Parameters used: Nk = 33, Ng = 512, Nh = 128 and Np = Ns = 64. 



B. 2D square lattice 

The bandwidth of the 2D square lattice is 8 and we have chosen U = 4 and [/ = 18 as prototypes for a 
metal and an insulator respectively. As in ID A^^ = 4 is sufficient. We will first discuss the metal case. The 
reasoning and organization follows closely the setup in the previous section. We have chosen the parameters 
V — 1.5 and V' = 0.75 in the metal case U ~ 4. In Figs. (|20I21|) the lattice Green's function arc shown 
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FIG. 20: Imaginary part of the site-diagonal Green's function at the F-point. Parameters used for DMFT-I-GW: 
Nk = 169, Ng = 512, Nh = 128 and Np ^ Ns = 64. 
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FIG. 21: Real part of the site-nondiagonal Green's function at the F-point. Parameters used for DMFT+GW: 
Nk = 169, Ng = 512, Nh = 128 and Np = Ns = 64. 



whereas the corresponding polarization and screened interaction are displayed in Figs. (|22I24() . At self- 
consistency the impurity screened interaction is 0.7 and the effective Hubbard U = 2.4, strongly reduced 
compared to the bare values. We stress that the amount of screening that are taking place is in general 
dependent of the non- locality parameters V and V' , which in this work is chosen arbitrarily. In 2D the 
large iuj limit is numerically satisfied: the diagonal part approaches AV' and the non-diagonal part AV 
respectively. It is worth mention that the overall magnitude of the polarization function P'~^^ [iuj] decreases 
for increasing U. As a consequence, the overall magnitude of the correlated part of the screened interaction 
W^liuj) increases for increasing U. As a comparison to the ID case, the screened interaction on real axis at 
the F-point is shown in Fig. EH 
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FIG. 22: Real part of the polarization function at the F-point. The static impurity contribution is P\{iu)m = 0) 
-1.1. Parameters used for DMFT-FGW: Nk = 169, Ng = 512, Nh = 128 and Np = Ns = 64. 
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FIG. 23: Real part of the site-diagonal (correlated) screened interaction at the F-point. The bare Hubbard U has 
been subtracted. Parameters used for DMFT+GW: Nk = 169, Ng = 512, A^^ = 128 and Np = Ns = 64. 
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FIG. 24: Real part of the sito-nondiagonal screened interaction at the F-point. Parameters used for DMFT-I-GW: 
A^fe = 169, Ng = 512, Nh = 128 and Np = Ns = 64. 
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FIG. 25: Real and imaginary part of the site-diagonal screened interaction at the F-point. We used an artificial 
broadening of 0.5. Parameters used for DMFT+GW: Nk = 169, Ng = 512, Nh = 128 and Np = Ns = 64. 



As an illustration we display in Fig. H26|l the GW self-energy, which clearly exhibits Fermi-liquid charac- 
teristics, derived using Eq. ^| 
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FIG. 26: Real and imaginary part of the GW self-energy including the double counting term (site 1 and spin up) 
for a typical metal. We used an artificial broadening of 0.5. Parameters used for DMFT-f GW: Nk = 33, Ng = 512, 
Nh = 128 and Np = Ns = 64. 



The metallic LDOS and a typical quasiparticle spectral function are shown in Figs. H27I28(I . The downward 
shift of the quasiparticle position is consistent with the scenario observed in the ID case. 
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FIG. 27: Local density of states for site 1 and spin up. The chemical potential corresponds to energy zero. We 
introduced an artificial broadening of 0.25. Parameters used for DMFT+GW: Nk = 169, Ng = 512, Nh = 128 and 
Np =Ns = 64. 




FIG. 28: Spectral function at the P-point. The chemical potential corresponds to energy zero. The solid line 
corresponds to the DMFT+GW case and the dashed line to the DMFT case. We used an artificial broadening of 0.5. 
Parameters used: A^^ = 169, Ng = 512, Nh = 128 and Np ^ Ns = 64. 



Let us finally consider the strong insulator case U = 18. We have chosen the parameters V — 4 and V — 3 
which can be considered as a substantial ofF-site interaction, however there exists no large difference in the 
DMFT Green's function compared to the DMFT+GW one (see Figs. (|29I30|) '). The impurity screening is 
found to W{iuJm = 0) = 16.9 the effective Hubbard U — 17.1, implying a reduced bandgap. 
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FIG. 29: Imaginary part of the site-diagonal Green's function at the F-point. Parameters used for DMFT+GW: 
Nk = 169, Ng = 512, Nh = 128 and Np ^ Ns = 64. 
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FIG. 30: Real part of the site-nondiagonal Green's function at the F-point. Parameters used for DMFT-(-GW: 
Nk = 169, Ng = 512, Nh = 128 and Np = Ns = 64. 



The corresponding spectral function is shown in Fig. H31|l . We note that in the strong insulator case the 
Hubbard gap U) is indeed somewhat reduced due to the inclusion of the GW self-energy. 
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FIG. 31: Spectral function at the F-point. The chemical potential corresponds to energy zero. The solid line 
corresponds to the DMFT+GW case and the dashed line to the DMFT case. We used an artificial broadening of 0.5. 
Parameters used: Nk = 169, Ng = 512, Nh = 128 and Np = Ns = 64. 



In the present study a full self- consistency is performed, including the non-local GW self-energy, in the local 
single site DMFT approach and the applicability of the method is tested for a model system. Eventually at 
self-consistency the full self-energy and polarization operator are obtained, from which e.g the full screened 
interaction is accessible. Far from the metal-insulator transition the combination of the GW method and 
the single site DMFT is from a numerical point of view fast and stable, even when a simple linear mixing 
scheme is utilized. Changes with respect to DMFT are in some cases substantial, and are related to the 
long-rangeness of the GW kernel, specified by two hopping parameters. 

Next we will study the 2D metal-insulator transition as well as doping away from half-filling. 



IV. CONCLUDING REMARKS 
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^® Due to the deflntion ofW'' = W-U the contribution 

S?j>(q) = - ^^^^ ^'^^ ^ Gijajk; ivn){cos{kx - qx) + cos{ky - qy)) 



^ ^ GijCT(k; Wn)cos{)tx - qx)cos{ky - qy) 



n k 

is implicitely included in E°. 

Symmetry for G and E in the 2D case: diagonal elements GuaCk) = GuaOi + G) for all G = n{ni + n2,n\ — n2). 
Non-diagonal elements [i / j) Gija{\<i) — Gijcr{k + G) if (ni + 112) even and Gijai^) = —Gijo- (k-|- G) if (ni +712) 
odd. Further Gijo-(k) = Gija{TZ\i), where k lies in the IBZ, because /i(k) — h(TZk). 

The polarization operator obeys P{m) = ^ G{n)G{n + m) = ^ G* {—n)G* {—n — m) = G* {n)G* {n — m) = 
E,, G{n)G{n - m)]* = P{-m)*. The self-energy operator obeys"E(n) = G{n - m)W{m) = G*{-n + 

'm,)VV'*(— m) = 

E„ G*(-n - m)W*{m) = G{-n - m)W(m)]* = E(-n)*. 
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The Kramer-Kronig (KK) relations can be written as 

ReE(u;) = -i /do^'ini^ (52) 

TV J LJ — LU' 

ImE(u.^-^ / ..'ReE(.') 



(ui) = — dw' ■ 
^ J 



OJ — UJ' 



which are very important in connecting the real and imaginary part of a given complex function. With the 
knowledge of the self-energy on the whole real axis on can perform an analytical continuation to the imaginary 
axis 



E(.) = -i fd.''^ 

TT J Z — Ul' 



for any z in the upper half-plane. It is readily shown that ReE(ia;) = and ImE(icij) = — ImE(— iw) if ImE(a)) 
symmetric around = 0. Furthermore it cam be shown that the slopes of ReE(a;) and ImE(ia;) are equal for small 
Ul and ImE(iw) ~ 5{u!) if the spectral function has zero amplitude at the chemical potential a; = 0. 
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